/* 
 * File:   Pattern.cpp
 * Author: Marco Ambu
 * 
 * Created on September 23, 2010, 11:01 AM
 */

#include "pr/Pattern.h"
#include "util/Matrix.h"
#include <algorithm>
#include <cmath>
#include <stdexcept>

namespace pr {

Pattern::Pattern(std::size_t featureCount)
  : pattern_(featureCount, 0.0)
{
}

Pattern::Pattern(const std::vector<util::Real>& values)
  : pattern_(values)
{
}

std::size_t Pattern::featureCount() const
{
  return pattern_.size();
}

util::Real Pattern::getFeature(std::size_t i) const
{
  return pattern_[i];
}

void Pattern::setFeature(std::size_t i, util::Real value)
{
  pattern_[i] = value;
}

util::Real Pattern::operator[](std::size_t i) const
{
  return pattern_[i];
}

util::Real& Pattern::operator[](std::size_t i)
{
  return pattern_[i];
}

util::Real Pattern::min() const
{
  return *std::min_element(pattern_.begin(), pattern_.end());
}

util::Real Pattern::max() const
{
  return *std::max_element(pattern_.begin(), pattern_.end());
}

void Pattern::maxpos(util::Real& max, std::size_t& pos) const
{
  const_iterator it = std::max_element(pattern_.begin(), pattern_.end());
  max = *it;
  pos = it - pattern_.begin();
}

util::Real Pattern::mean() const
{
  util::Real m = 0.0;

  for (const_iterator i = pattern_.begin(); i != pattern_.end(); ++i)
    m += *i;
  return m/pattern_.size();
}

util::Real Pattern::variance() const
{
  util::Real var = 0.0, m = mean();

  for (const_iterator i = pattern_.begin(); i != pattern_.end(); ++i)
    var += (*i - m)*(*i - m);
  return var/(pattern_.size()-1);
}

void Pattern::meanAndVariance(util::Real& meanValue, util::Real& variance)
{
  meanValue = mean();

  variance = 0.0;
  for (const_iterator i = pattern_.begin(); i != pattern_.end(); ++i)
  {
    double diffFromMean = *i - meanValue;
    variance += diffFromMean * diffFromMean;
  }
  variance /= (pattern_.size() - 1);
}

util::Real Pattern::median() const
{
  std::vector<util::Real> copy(pattern_.begin(), pattern_.end());
  std::sort(copy.begin(), copy.end());
  return copy[pattern_.size()/2];
}

util::Real Pattern::normInfinite() const
{
  return *std::max_element(pattern_.begin(), pattern_.end());
}

util::Real Pattern::norm1() const
{
  util::Real n = 0.0;

  for (const_iterator i = pattern_.begin(); i != pattern_.end(); ++i)
    n += fabs(*i);
  return n;
}

util::Real Pattern::norm2() const
{
  util::Real n = 0.0;

  for (const_iterator i = pattern_.begin(); i != pattern_.end(); ++i)
    n += (*i) * (*i);
  return sqrt(n);
}

Pattern& Pattern::normalizeInfinite()
{
  util::Real normInverse = 1.0 / normInfinite();

  for (iterator i = pattern_.begin(); i != pattern_.end(); ++i)
    *i *= normInverse;  // multiply by 1/norm
  return *this;
}

Pattern& Pattern::normalize2()
{
  util::Real n = norm2();
  if (!(n > 0.0 || n < 0.0))
    return *this;

  n = 1.0 / n;
  for (iterator i = pattern_.begin(); i != pattern_.end(); ++i)
    *i *= n;  // multiply by 1/norm
  return *this;
}

Pattern& Pattern::normalizeByMaxRange()
{
  util::Real minimum = min();
  util::Real range = max() - minimum;

  for (iterator i = pattern_.begin(); i != pattern_.end(); ++i)
    *i = (*i-minimum) / range;
  return *this;
}

Pattern& Pattern::normalizeGauss()
{
  util::Real m, sigma2;
  meanAndVariance(m, sigma2);

  for (iterator i = pattern_.begin(); i != pattern_.end(); ++i)
    *i = (*i-m) / sigma2;
  return (*this);
}

Pattern& Pattern::sort(bool reverse)
{
  if (reverse)
    std::sort(pattern_.rbegin(), pattern_.rend());
  else
    std::sort(pattern_.begin(), pattern_.end());
  return *this;
}

Pattern& Pattern::operator+=(const Pattern& pattern)
{
  if (featureCount() != pattern.featureCount())
    throw std::invalid_argument("Cannot sum patterns with different size");

  const_iterator it = pattern.pattern_.begin();
  for (iterator i = pattern_.begin(); i != pattern_.end(); ++i, ++it)
    *i += *it;
  return *this;
}

Pattern& Pattern::operator-=(const Pattern& pattern)
{
  if (pattern_.size() != pattern.pattern_.size())
    throw std::invalid_argument("Cannot subtract patterns with different size");

  const_iterator it = pattern.pattern_.begin();
  for (iterator i = pattern_.begin(); i != pattern_.end(); ++i, ++it)
    *i -= *it;
  return *this;
}

Pattern& Pattern::operator+=(util::Real value)
{
  for (iterator i = pattern_.begin(); i != pattern_.end(); ++i)
    *i += value;
  return *this;
}

Pattern& Pattern::operator-=(util::Real value)
{
  for (iterator i = pattern_.begin(); i != pattern_.end(); ++i)
    *i -= value;
  return *this;
}

Pattern& Pattern::operator*=(util::Real value)
{
  for (iterator i = pattern_.begin(); i != pattern_.end(); ++i)
    *i *= value;
  return *this;
}

Pattern& Pattern::operator/=(util::Real value)
{
  if (!(value > 0.0) && !(value < 0.0))
    throw std::invalid_argument("Error: division by 0");

  util::Real valueInverse = 1.0 / value;
  for (iterator i = pattern_.begin(); i != pattern_.end(); ++i)
    *i *= valueInverse;
  return *this;
}

Pattern& Pattern::operator^=(util::Real value)
{
  for (iterator i = pattern_.begin(); i != pattern_.end(); ++i)
    *i = pow(*i, value);
  return *this;
}

util::Real Pattern::euclideanDistance(const Pattern& pattern) const
{
  if (featureCount() != pattern.featureCount())
    throw std::invalid_argument("Cannot compute distance between to patterns with different size.");

  util::Real dist = 0.0;
  const_iterator i = pattern_.begin(), pi = pattern.pattern_.begin();
  for ( ; i != pattern_.end(); ++i, ++pi)
    dist += (*i - *pi)*(*i - *pi);

  return sqrt(dist);
}

util::Real Pattern::weightedEuclideanDistance(const Pattern& pattern, const std::vector<util::Real>& weight) const
{
  if (featureCount() != pattern.featureCount())
    throw std::invalid_argument("Cannot compute distance between to patterns with different size.");

  util::Real dist = 0.0;
  const_iterator i = pattern_.begin(),
          pi = pattern.pattern_.begin();
  std::vector<util::Real>::const_iterator wi = weight.begin();
  for ( ; i != pattern_.end(); ++i, ++pi, ++wi)
    dist += (*i - *pi)*(*i - *pi)*(*wi);

  return sqrt(dist);
}

util::Real Pattern::cityblockDistance(const Pattern& pattern) const
{
  if (featureCount() != pattern.featureCount())
    throw std::invalid_argument("Cannot compute distance between to patterns with different size.");

  util::Real dist = 0.0;
  const_iterator i = pattern_.begin(),
          pi = pattern.pattern_.begin();
  for ( ; i != pattern_.end(); ++i, ++pi)
    dist += fabs(*i - *pi);

  return dist;
}

util::Real Pattern::mahalanobisDistance(const Pattern& pattern, const std::vector<util::Real>& weight) const
{
  if (featureCount() != pattern.featureCount())
    throw std::invalid_argument("Cannot compute distance between to patterns with different size.");

  util::Real dist = 0.0;
  const_iterator i = pattern_.begin(),
          pi = pattern.pattern_.begin();
  std::vector<util::Real>::const_iterator wi = weight.begin();
  for ( ; i != pattern_.end(); ++i, ++pi, ++wi)
    dist += (*i - *pi)*(*i - *pi)/(*wi);
//    dist -= (*i)*(*pi)/sqrt(*wi); // -= !!!

  return sqrt(dist);
}

 util::Real Pattern::scalarProduct(const Pattern& pattern) const
 {
  if (featureCount() != pattern.featureCount())
    throw std::invalid_argument("Cannot compute scalar product between to patterns with different size.");

  util::Real  prod = 0.0;
  const_iterator i = pattern_.begin(),
          pi = pattern.pattern_.begin();
  for ( ; i != pattern_.end(); ++i, ++pi)
    prod += (*i * *pi);

  return prod;
}

 void Pattern::selfProduct(util::Matrix& product)
 {
   if (product.rows() != featureCount() || product.cols() != featureCount())
     throw std::invalid_argument("Output matrix not valid for self product.");

   for (std::size_t i = 0; i < pattern_.size(); ++i)
    for (std::size_t j = 0; j < pattern_.size(); ++j)
      product(i, j) = (*this)[i]*(*this)[j];
 }

// ATTENZIONE, questa proiezione modifica la dimensione del vettore
Pattern& Pattern::project(const std::vector<Pattern>& vp)
{
  *this = projection(vp);
  return *this;
}

Pattern Pattern::projection(const std::vector<Pattern>& vp) const
{
  Pattern proj(vp.size());
  unsigned int i;

  for (i=0; i < vp.size(); ++i)
    proj[i] = scalarProduct(vp[i]);
  return proj;
}


Pattern operator+(const Pattern& left, const Pattern& right)
{
  Pattern result(left);
  result += right;
  return result;
}

Pattern operator-(const Pattern& left, const Pattern& right)
{
  Pattern result(left);
  result -= right;
  return result;
}

Pattern operator+(const Pattern& left, util::Real value)
{
  Pattern result(left);
  result += value;
  return result;
}

Pattern operator-(const Pattern& left, util::Real value)
{
  Pattern result(left);
  result -= value;
  return result;
}

Pattern operator*(const Pattern& left, util::Real value)
{
  Pattern result(left);
  result *= value;
  return result;
}

Pattern operator/(const Pattern& left, util::Real value)
{
  if (!(value > 0.0) && !(value < 0.0))
    throw std::invalid_argument("Error: division by 0");

  Pattern result(left);
  result /= value;
  return result;
}

Pattern operator^(const Pattern& left, util::Real value)
{
  Pattern result(left);
  result ^= value;
  return result;
}

} // namespace pr
